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ABSTRACT 

We constrain the evolution of the rest-frame far-infrared (FIR) luminosity function 
out to high redshift, by combining several pieces of complementary information pro- 
vided by the deep Balloon-borne Large- Aperture Submillimeter Telescope surveys at 
250, 350 and 500 |im, as well as other FIR and millimetre data. Unlike most other 
phenomenological models, we characterise the uncertainties in our fitted parameters 
using Monte Carlo Markov Chains. We use a bivariate local luminosity function that 
depends only on FIR luminosity and 60-to-100 |imi colour, along with a single library 
of galaxy spectral energy distributions indexed by colour, and apply simple luminosity 
and density evolution. We use the surface density of sources, Cosmic Infrared Back- 
ground (CIB) measurements and redshift distributions of bright sources, for which 
identifications have been made, to constrain this model. The precise evolution of the 
FIR luminosity function across this crucial range has eluded studies at longer wave- 
lengths (e.g., using SCUBA and MAMBO) and at shorter wavelengths (e.g., with 
Spitzer), and should provide a key piece of information required for the study of 
galaxy evolution. Our adoption of Monte Carlo methods enables us not only to find 
the best-fit evolution model, but also to explore correlations between the fitted param- 
eters. Our model-fitting approach allows us to focus on sources of tension coming from 
the combination of data- sets. We specifically find that our choice of parameterisation 
has difficulty fitting the combination of CIB measurements and redshift distribution 
of sources near 1 mm. Existing and future data sets will be able to dramatically im- 
prove the fits, as well as break strong degeneracies among the models. Two particular 
examples that we find to be crucial are: obtaining robust information on redshift dis- 
tributions; and placing tighter constraints on the range of spectral shapes for low 
luminosity (Lfir < 1O 1O L ) sources. 

Key words: galaxies: evolution - galaxies: high-redshift - stars: formation - submil- 
limetre: galaxies 



1 INTRODUCTION 

It is now known that a significant fraction of the total light 
produced by stars and active galactic nuclei (AGN) through- 
out cosmic history is absorbed by dust and re-radiated ther- 
mally at much longer wavelengths. This light was first ob- 
served at low angular resolution by the COBE satellite as 
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the diffuse Cosmic Infrared Background (CIB, Puget et al 
1996; Fi xsen et al.||1998 ). Over the last decade it has been 



largely resolved into point sources at wavelengths ~ 100- 
1000 um, demonstrating that it is predominantly produced 



by individual galaxies ( 


Dole et al. 


2006 


Pope 2007 Serjeant 


et al.||2008 Marsden et al.||2009 


Pascale et al.||2009). Sur- 



veys with the IRAS, ISO and Spitzer satellites show that 
most of the shorter-wavelength light is produced by galax- 
ies at redshifts z < 1, while ground-based submillimetre 
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(submm) surveys have found most of the longer-wavelength 
light to be produced by more distant objects. Recent sur- 
veys by the Balloon-borne Large Aperture Submillimeter 
Telescope (BLAST) at 250, 350 and 500 ^im, a precursor to 
the Herschel/ SPIRE surveys that are now well underway, 
have shown that the transition from low to high redshifts as 
one observes at longer wavelengths occurs gradually across 
the 250-500 |xm band. 

The fact that there is a transition to higher-redshift 
sources observed at longer wavelengths is not a surprise. 
Many groups have predicted this general behaviour using 
simple parameterised models for the evolution of local far- 
infrared (FIR) and submm galaxy luminosity functions. The 
data typically fit at these wavelengths include the surface 
density of sources as a function of brightness (source counts) 



and redshift information, when available (e.g. 


, Blain & Lon- 


gair|1993||Guiderdoni et al.|1997||Blain et al. 


1999| ICharyfc 


Elbaz|2001||Malkan & Stecker|2001| |Rowan-F 


lobinson|2001| 


Lagache et al.|2003||Dole et al.|2003||Lagache et al.|2004||Le| 


Borgne et al. 2009; Valiante et al.|2009|). These phenomeno- 



logical models may be thought of as the simplest fitting func- 
tions available, since they typically include only two main 
ingredients: spectral energy distribution (SED) templates, 
to relate observed flux densities in different bands given lu- 
minosities and redshifts; and some evolutionary form for the 
luminosity function, to produce greater numbers of objects 
at higher redshifts - typically luminosity or density evolu- 
tion following a power law in (1 + z). 

As observational data at ~ 10-1000 \im have improved, 
in terms of wavelength coverage, survey area and depth, 
many authors have added greater complexity to their mod- 
els. For example, it is now common to divide up the local 
luminosity function into multiple galaxy populations, assign- 
ing different SEDs to each, and then evolving the popula- 
tions independently (e.g., |Rowan- Robinson] |2001 1 |Lagache] 



et al.||2003 ), or assuming some relation between the IR lu- 
minosity and the AGN content (Vali ante et al.|2009 ). With 
this added freedom, such models can simultaneously fit the 
longer-wavelength 850-1200 \im SCUBA/MAMBO/AzTEC 
number counts and approximate redshift distributions, as 
well as the shallower IRAS/ ISO /Spitzer surveys in the FIR 
(60-200 |^m), and more recently the deep 24 |xm Spitzer sur- 
veys (e.g.,|Lagache et al.|2004| | Valiante et al. 2009; Rowan 



Robinson|2009|). Despite these successes, we note that, prior 



to the first measurements, predictions for the number counts 



in the BLAST/SPIRE bands varied widely (Patanchon et al. 
2009 hereafter P09). We believe there are two main reasons 



for these discrepancies. First, the number of parameters as- 
sociated with the multiple discrete populations is large, po- 
tentially leading to significant uncertainties in any part of 
the spectrum that lacks observational constraints. Second, 
it is common knowledge that rest-frame dust temperature, 
and hence bolometric luminosity and total dust mass, are 
degenerate with source redshift (e.g., |Blain et al. 2003); a 
'redder' object could either be a cooler local galaxy, or a 
warmer, more distant galaxy. For this reason, assumptions 
about the SED shapes for each galaxy population, and the 
potential for evolution in these shapes, can significantly af- 
fect the results. 

An alternative to phenomenological modelling of the 
data is the ab initio approach, or solutions to the forward 
problem: simulate as much of the physics of galaxy forma- 



tion and evolution as possible, and evolve the model forward 
in time from the Big Bang, tweaking model parameters to fit 
observables (e.g., |Baugh et al. 2005). Currently, such mod- 
els usually incorporate the numerical evolution of the dark 
matter distribution, and adopt a range of recipes to assign 
galaxies to over-densities in this distribution, both as a func- 
tion of mass and the dark matter halo merger histories. Ob- 
servational constraints include the cosmic microwave and 
infrared backgrounds, the luminosity functions and spatial 
clustering of local galaxy populations, and more recently, 
their observed surface density and redshift distributions. 

For the purpose of understanding galaxy evolution the- 
ory, ab initio models presently offer the most complete tool- 
box, combining a wide variety of physics into a single coher- 
ent model. However, they are still quite limited in the preci- 
sion of inference that they allow, since the simulations lack 
the necessary resolution, and therefore assumptions about 
smaller-scale physics must be made (such as the details 
of star formation within molecular clouds, growth of black 
holes, and the interaction between the galaxy and the inter- 
galactic medium). These assumptions are guided by intu- 
ition and simple tuning recipes based, for example, on in- 
formation from higher-resolution hydrodynamic simulations 
(e.g., Narayanan et aL||2010 ). These limitations result in a 
large number of parameters (undoubtedly with many par- 
tial degeneracies which can be tuned to match the sub-grid 
physics) , and it can be extremely difficult to estimate mean- 
ingful uncertainties. 

In this paper we consider a phenomenological model 
with more modest goals. Unlike ab initio models, and more 
recent multi-population phenomenological models that seek 
to fit the widest range of data available (e.g., attempting to 
connect the submm to FIR and mid-IR galaxy populations, 
as in JChary fe Elb az 2001; |Lagache et "ah] [20041 |Valiante| 
et al.| 2009), we focus our efforts on data that constrain only 
the evolution of the rest-frame FIR peak (A ~ 60-200 urn), 
for galaxies at redshifts up to z ~ 4.5. This redshift range 
encompasses the bulk of the 850-1200 \im submm galaxy 
(SMG) population which peaks near z ^ 2.5 ( Aretxaga et al 



2003| |Chapman et al.||2003| |2005| |Chapin et al.||2009| ), as 



well as the most distant spectroscopically-connrmed exam- 
pies jCapak et al.|2008||Schinnerer et al.|2008||Daddi et al. 



2009 Coppin et al. 2009). We are therefore only attempt- 
ing to fit data that directly probe dust-reprocessed radiation 
from the most active star- forming galaxies, from their for- 
mation epoch to the present day. We also explicitly set out to 
determine whether a single galaxy population with a simple 
evolutionary form can reproduce the observed data across 
the peak in the CIB. Most authors have concluded, through 
fitting 'by hand', that multiple populations with indepen- 
dent evolutionary forms are required by the data; however, 
an exhaustive non-linear search of parameter space has never 
been performed as in this work to determine: (i) whether a 
more complicated model is indeed necessary; and if so (ii) 
identify precisely where the tension is coming from to probe 
the types of new models, and/or new data that would be 
required to fit such models. 

For local galaxies, the region of the spectrum we are 
considering is quite smooth. At wavelengths J> 100 \im, SED 
shapes resemble modified blackbodies, 
although at slightly shorter wavelengths they are typically 
brighter and flatter due to a combination of opacity effects, 
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and sensitivity to warmer dust grains. We therefore do not 
need to pay special attention to tuning the MIR spectra of 
our models, e.g., stochastically heated small dust grains, in- 
cluding polycyclic aromatic hydrocarbons, or PAHs (Draine 



&; Li 2001), as has been necessary to fit the deep Spitzer 
24 urn data flChary fe Elbaz||2001| |Lagache et al.|[2004 



Valiante et al. 1120091 |Rowan-Robinson||2009| . We do, how- 



ever, include a range of SED shapes in our analysis, with a 
distribution characterised by the single C = log 10 (5 f 6o/*S'ioo) 
colour near the peak of the rest-frame FIR emission (which 
is a good indicator of the FIR peak wavelength) and its 
well-known correlation with FIR luminosity, L = Lfir, the 
integrated luminosity from 42.5-1 22.5 urn, (|Soifer &; Neuge-| 
|bauer||199lT|Chapman et al.|2003| |Chapin et al.|2009| ). 



We combine this simple SED shape parameterisation 
with the local IRAS luminosity function as our local bound- 
ary condition. We then evolve this local bivariate distribu- 
tion, $(L,C), using only luminosity and density evolution, 
to fit the submm-FIR data. While the idea of incorporat- 
ing a correlation between luminosity and FIR colour is not 
new (e.g., Lagache et al. 2003| Lewis et al.||2005] Valiante 



et al. 2009), our reduced number of parameters allows us 



to fully explore the parameter space using Monte Carlo 
Markov Chains (MCMC). To our knowledge, there are only 
two other published attempts to fully characterise the uncer- 
tainties in a phenomenological model - Kell y et al.| ([2008 ) 



and Le Borgne et al. ( 2009| , only the latter of which was 



concerned with submm-FIR surveys. In addition, a concur- 



rent study by Bethermin et al. (2011) uses methods similar 



to those presented here. Another feature of our analysis that 
sets it apart from earlier work is our focus on the potential 
evolution in the correlation between luminosity and colour, 
since it is degenerate with redshift and heavily influences 
conclusions about the dust-enshrouded star-formation rate 
history. 

Throughout this paper, we consider two basic models: 
one in which the local correlation holds at high redshift; 
and a second in which the correlation undergoes luminosity 
evolution (making galaxies of a given luminosity in the past 
appear cooler than at the present day). This is an area for 
which BLAST data, and newer SPIRE data, provide the 
strongest constraints on this crucial part of the spectrum. 
This approach has allowed us both to: choose the model 
that best-fits existing data; and clearly indicate what future 
data are required to break the remaining degeneracies in 
parameter space. 

A basic assumption that we make, as with all mod- 
els of this type, is that high-redshift luminosity functions 
smoothly evolve over time to produce the modern-day z = 
luminosity functions. If there is a significant galaxy popu- 
lation that existed in the early Universe, but is completely 
absent in the local Universe (even as a faint tail), our model 
will not give plausible results. We also emphasise the fact 
that our model will not give useful predictions for data far 
from the rest-frame FIR peak, such as the 24 urn source 
counts. Further work would be needed in order to achieve 
this, and in all likelihood, require additional model param- 
eters. 

In Sections |2. 1 1 and [272] we describe the local boundary 
conditions of our model - the local luminosity and colour- 
luminosity distributions, and our adopted SED templates. 
The parameterisation of the redshift evolution is described 



in Section [2. 3| and the connection of this model to observed 
quantities (such as number counts, background intensities 
and redshift distributions) is provided in Section |2.4| The 
data sets that we use to fit the model are given in Section [3] 
and our fitting procedure is summarised in Section [4] The 
results of the fits are presented in Section [5] and Section [6] 
discusses the implications of and future improvements to the 
model. 

Throughout this paper a standard cosmology is 
adopted, with Q M 0.272, fi A .728, and H 
70.4kms _1 Mpc -1 (Komatsu et al.|2011 ). 



2 EVOLUTION OF THE LUMINOSITY 
FUNCTION 

2.1 Local Luminosity and Colour Distributions 

In most past phenomenological models, authors have used 
either the local 60 urn luminosity function (e.g., Saunders 
et al.|1990|) or the lo cal 850 um SCUBA luminosity function 



flDunne et aT]|2000| ). The former is one of the most well- 
studied luminosity functions, based on the all-sky IRAS 12, 
25, 60 and 100 um survey. Since SCUBA was not sensitive 
enough to conduct a survey over a significant portion of the 
sky, pointed follow-up of an IRAS galaxy sample was em- 
ployed for the latter. This technique is adequate, provided 
that no significant population of cool (< 25 K) galaxies ex- 
ist in the local Universe. BLAST traces the peak of dust 
emission at high redshift, as IRAS does locally; emission 
at 850 um comes from the Ray leigh- Jeans part of the spec- 
trum. For these reasons, we choose an IRAS-based luminos- 
ity function as the basis for our model. 

Since we are interested in the range of SED shapes 
that produce both the rest-frame submm and FIR emis- 
sion, we also use the distribution of C = log^Seo/Sioo) 
colours as a function of luminosity. The observed correlation, 
with no corrections for observational biases, was measured 



by Soifer &; Neugebauer (1991) and has been used in some 
phenomenological models (e.g., Lagache et al. 120031). An at- 



tempt was made to measure the full bivariate luminosity- 
colour distribution using the 1/V r 
|1968| > by |Chapman et al.| p003[ 



max technique (Schmidt 
However, we instead use 



the updated version of this distribution from Ch apin et al.| 
(2009, hereafter C09a), which incorporates additional cor- 
rections for the IRAS bandpasses, a bias against detecting 
cooler galaxies in the original 60 um flux-limited sample, and 
redshift evolution. We make one minor alteration to the dis- 
tribution; as noted in C09a (and other previous authors, 
e.g. Saunders et al. 1990; Lawrence et al.p999| , the faint- 
end of the luminosity function is biased high by the local 
over-density of galaxies. Since the joint density of galax- 
ies as a function of luminosity and colour is formulated in 
C09a as <E>(L,C) = $(L)p(C\L) (i.e., the product of a pure 
luminosity function, with the conditional probability of a 
galaxy having a colour C given a luminosity L), we simply 



replace &(L) with the measurement from | Saunders et al 
( |1990| ), which is valid because they used an estimator that 
is insensitive to this over-density (and which was shown in 
C09a to be consistent at luminosities L > L*). Here (and 
throughout), L is defined to be the integrated 42.5-122.5 um 
FIR luminosity. Throughout, we follow the convention that 
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<E>(L) is the number density of objects per unit luminosity 
and 4>(L) is the number density per decade of luminosity. 



2.2 SED Library 

As mentioned in the introduction, the shapes of the submm- 
FIR SEDs of most galaxies in the local Universe are 
reasonably-well parameterised by a simple 2-parameter 
modified blackbody (in addition to the normalisation), S u oc 
i/ p B(i/,T), where B(y,T) is the Planck function. While a 
range of values of /3 have been measured, they typically 
do not vary much from a canonical value of f3 = 1.5. Fur- 
thermore, /3 and T G bs are highly anti-correlated in the fits. 
Therefore, it is plausible that a single parameter can accu- 
rately describe most of the spread in the shapes of locally- 
observed FIR SEDs - in our case we use the FIR colour C. 
In further support of this simple parameterisation, |Dunne| 
et al. ( 2000 ) showed that it was possible to map the 60 um 



luminosity function to the 850 um luminosity function by 
adopting SEDs that follow the observed correlation between 
temperature and luminosity. Similarly, Serjea nt &; Harrison| 
(2005) fit temperatures to IRAS 60 and 100 um data with 
a fixed /3 — 1.3, and for each object estimated their 850 um 
brightnesses. They found, using these predicted 850 um flux 
densities, that they could also map the IRAS luminosity 
function to the 850 um luminosity function. 

In the spirit of these earlier analyses, we seek a set 
of SED templates that can transform the IRAS colour- 
luminosity distribution to local luminosity functions at other 



adjacent wavelengths, namely: IRAS 12 urn (Fang et al. 
; ISOCAM 15 um ( |Xu|2000|); IRAS 25 um Jshupe et al. 



1998 



1998 



; and SCUBA 850 urn ( [Dunne & Eales||2Q0l| ) | x | The 
shortest wavelength data that we will attempt to fit are 
the observed source counts at 70 um. It is therefore im- 
portant to match the local luminosity functions at wave- 
lengths as short as 12 urn, since this wavelength is redshifted 
to 70 um at z ~ 5 (we only expect a minor contribution to 
the 12-15 um SEDs from more complicated emission mecha- 
nisms, and this will only impact the 70 um measurements for 
z ^ 4.5 sources). On the long- wavelength side, by achieving 
consistency between the IRAS and 850 um luminosity func- 
tions, we can expect to reasonably interpolate the rest-frame 
luminosity functions in the BLAST bands that are brack- 
eted by these wavelengths. Again, we emphasise that we are 
not attempting to fit 24 um number counts, which would re- 
quire accurate modelling of the data at shorter wavelengths 
(<^ 8 um) where the SEDs are considerably more compli- 
cated. 

We have examined and rejected three common SED 
models for individual galaxies. First, we attempted to use 
single-temperature modified blackbodies. We produced a li- 
brary of SEDs with fixed values of /3 and a range of tem- 
peratures. We then predicted each of the monochromatic 
luminosity functions, <f>(L u ), from <J>(L,C), 



(L(L V ,C),C) 



dL(L u ,C) 



dC, 



(1) 



1 Since each study adopted slightly different values of Ho, we 
have corrected the luminosities and volumes for the value Ho = 
70.4 km s — 1 Mpc - 1 used in this paper. 



where L(L U ,C) is the FIR luminosity for an SED in our 
library of colour C, normalised to the luminosity density 



L u at a frequency v. Similar to Dunne et al. (2000) and 



Serjeant h Harrison (2005), we were able to obtain good 



agreement between the IRAS and 850 um luminosity func- 
tions using /3 = 1.5. However, unsurprisingly, this simplistic 
model is a poor fit to the shorter wavelength data, since 
the Wien tail of our single-temperature SEDs falls off con- 
siderably more rapidly than for real galaxies, which con- 
tain a mixture of dust at different temperatures and com- 
positions. Even the apparent plausibility of our SEDs for 
the longer- wavelength data can be deceiving; an effective 
P — 1.5 single-temperature spectrum can also be produced 
by the superposition of a range of dust populations with a 
steeper f3 = 2.0 at a range of appropriately-selected temper- 
atures. This fact serves to remind us that no simple physical 
meaning should be attached to these model parameters; the 
modified blackbody is only a convenient fitting function. 

Next, we tested two more realistic SED libraries (each 
spanning the submm-IR wavelengths of interest) that are 
commonly used in the submm and FIR literature: the tem- 



plates from Chary Elbaz| (2001) that were fit to data 



spanning 0.44-850 um, and which were used in their phe- 
nomenological model, based on evolution of the 15 um lumi- 
nosity function; and the templates of star-forming galaxies 
fro m|Dale et al.|(|2001| tha t were fit to IRAS and ISO data. 
The |Chary Elbaz| ( |20Ql| ) templates provided a reasonable 
extrapolation to the luminosity functions at 12-25 um, but 
led to a significant over-prediction at 850 um. The Dale et al.| 
(2001 J SEDs performed much better at 850 um, but led to 
moderate over-predictions at 12-25 um. Due to these short- 
comings, we decided to produce our own SED templates that 
vary smoothly as a function of C. 

The basis of our SED library is the model of |Draine| 



& Li (2007). Their parameterisation includes: (i) a set of 
mid-IR templates as a function of the PAH abundance, 
qpAH'i and (ii) longer- wavelength templates for cooler ther- 
mal emission that is composed of dust heated both by a 
single low-intensity radiation field, U m in, and a second com- 
ponent heated by a range of radiation intensities from Umin 
to [/max, where a factor 7 gives the fraction of the total 
dust emission produced by this second component. We ex- 
perimented with these 4 parameters (#pah, Umin, Umax and 
7) to produce a sequence of smoothly varying SEDs as a 
function of C that resulted in good extrapolations to the 
monochromatic luminosity functions on both sides of the 
FIR peak. 

We found that the results did not depend particularly 
heavily on (/pah and so thereafter we simply fix it to an 
intermediate value of 2.50 (from a possible range spanning 
0.10-4.58). The value of Umin effectively sets the apparent 
temperature of the coolest dust, and larger values of Umax 
and 7 increase the temperature and fraction of the hotter 
dust (i.e., together these parameters control most of the 
submm-FIR SED shape). We obtained good fits for our ex- 
trapolated luminosity functions by fixing Umax = 10 4 (from 
a possible range spanning 10 3 -10 6 ), and stepping through 
all 22 of the supplied templates, corresponding to values of 
0.1 ^ Umin ^ 25.0. We simultaneously increased 7 logarith- 
mically from 10 -9 to 0.4 over the 22 levels. Finally, for this 
set of SEDs, we found that the values of C only spanned 
—0.50 to 0.10, thereby missing some of the warmest values 
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Figure 1. SEDs used in our model, generated from a range of 
templates presented in |Draine fe Li| (2007). The 60-to-100 urn 
colour C of these SEDs ranges from —0.5 to 0.5, indicated by 
a colour gradient from red (cool, C = —0.5) to blue (warm, C = 
0.5). The abrupt increase in the density of SEDs which peak at 
A ^ 30 um is caused by an extrapolation from the warmest |Draine| 
[FTi] (2007) model (with C = 0.1) to even larger values of C 
by adding modified blackbodies with (3 = 1.5 and temperatures 
ranging from 47-98 K. A radio component is added on, based on 
the FIR-radio correlations; this dominates at A > 0.5-2 mm. 
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Figure 2. The local luminosity density functions derived from 
our adopted local luminosity function <£(£), colour-luminosity 
function p(C\L) and SED library, described in Secti ons |2.1| and 
[2.2| The 12, 15, 25 and 850 urn data sets are from |Fang et al.| 
(1998t, pCuH2005t , |Shi^>e et al.| ( fl998| ) and punne et al.| ( |2000| > 
respectively. 



qm log 



Sir/3.75 x 10 12 
Wm- 2 



log ( 



Wm- 2 Hz 



(2) 



and also a power law index a ra d- We add the radio con- 
tinuum using giR = 2.41 and a ra d = —0.4 at wavelengths 
longer than 100 p,m|^]This has a small effect on the number 
counts predictions at millimetre wavelengths, presented in 
Section [5] 

Fig. [2] shows the derived local monochromatic luminos- 
ity functions using the above SED library. The good cor- 
respondence indicates that we have found a plausible 2- 
parameter model for the distribution of galaxy luminosities 
and colours spanning 12-850 ^im. Later, in Section |6J"} we of- 
fer further evidence that this SED library is reasonable over 
the wavelength range of interest by comparing the colours of 
galaxies sampled from our best-fit model with real surveys 
that span 24-850 |^m. 



identified in the local Universe, and to a lesser extent, some 
of the cooler values (see fig. 4 in C09a). Since the Rayleigh- 
Jeans side of the Draine &; Li| ([2007 ) models resemble mod- 
ified blackbodies with /3 = 1.5, we simply extended the li- 
brary to a larger value of C = 0.5 by taking the C = 0.10 
SED and adding on modified blackbodies with ft — 1.5, tem- 
peratures T ranging from 47-98 K, and normalised such that 
they pass through the C = 0.10 model at 100 \im. At the red 
end, we simply assign the C = —0.50 SED to galaxies with 
C < —0.50 (i.e., galaxies that appear to have temperatures 
T <J 27 K using a f3 = 1.5 modified blackbody are truncated 
at T = 27 K). Fig. [I] shows the complete set of SEDs. The 
transition to the extrapolated C > 0.10 SEDs are obvious 
as an abrupt increase in the density of templates that peak 
at wavelengths A <J 30 \im. 

While it does not affect the results of the model fits, 
we also add on a power- law radio component to the SEDs, 



based on the FIR-radio correlation of Ivison et al. ( 2010 ) 



These authors measure the quantity qm, the ratio of the 
rest-frame 8-1000 \im flux to the 1400 MHz flux density, 



2.3 Extension to High Redshift 

We use a simple extension of the local luminosity function 
$(L, C) to high redshift, incorporating parametric forms for 
density p(z), luminosity g(z) and colour-luminosity h(z) evo- 
lution: 



&(L,C,z) = p{z) x $ 



9(z) 



xplC 



L 

W) 



(3) 



Here, <£(£, C, z) is the comoving luminosity evolution func- 
tion, with units Lq 1 Mpc -3 . Then the number of galaxies in 
a cell centred at (L, C, z) is: 



N = $(L, C, z) AL AC ^ Az, 

dz 



(4) 



2 The quantity qm has been updated to 2.40 using Her- 
schei/SPIRE observations pvison et al.|2010] >, but as the change 
is significantly smaller than the measurement errors, we have not 
updated our value. 
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where AL, AC and Az are the dimensions of the cell. 

We note that if &(L) were a power law, p and g would 
be completely degenerate (this is a well-known problem, 
e.g., Saunder s et al.| [l990). Our &(L) has both a break 
(at Li = 4 x 10 y L ) and a low- luminosity cut-off (at 
Lcut = 10 8 L©). These features, combined with our use 
of a colour- luminosity correlation, serve to break some of 
the degeneracy between the parametric functions, although 
strong correlations remain. We have not fully tested the ef- 
fect of v arying the low-luminosity cut-off (although see Sec- 
tion 6.5), but we note that the total luminosity, j L&(L)dL, 
integrated from L cut to infinity, is ~ 90 per cent of the full 
integral. 

To reduce the number of free parameters, we set 
h(z) = [g(z)] a . The parameter a controls the amount of 
colour- luminosity evolution in the model: with a = 1, the 
right-hand side of Equation [3] can be written as p(z) x 
&(L/g(z), C); with a = 0, the colour-luminosity relation 
does not evolve with redshift. The parameterisation of p(z) 
and g(z) is discussed in Section 4.1 We explore models with 



a fixed to 1.0 and 0.0, as well as with a allowed to vary as a 
free parameter. In all cases, p and g are free parameterised 
functions. The consequences of varying a are significant and 
are discussed in Section l(T4l 



2.4 Observables 

Our evolving luminosity function can be integrated across 
the appropriate variables to provide observables from the 
model. We first change variables from intrinsic luminosity L 
to observed flux density S u : 

/(&,, C,z) = * {L{S V , C, z), C, z) (H) (^) , (5) 
with 

L{S„,C,z) = - ■ 4 :^ Z ]f: , „ (6) 



(l + z)T(C,(l + z)vY 

where Dl is the luminosity distance and T(C, v) converts 
FIR luminosity L to luminosity density L v at rest-frame 
frequency v for the SED template with colour C. The dV/dz 
term converts the counts from number per unit volume to 
number per unit redshift. 

We use four types of data in this analysis: 

• differential number counts, calculated by integrating 
across colour and redshift. 



dN(Su) 
dSu 



f(S u ,C, z) dCdz- 



• background intensity (CIB), obtained by further inte- 
gration over S u , 



poo 
/ Su 

Jo 



dN(Su 



dSu 



dSu] 



(8) 



• background intensity as a function of redshift of emit- 
ting sources, 



dlu 

dz 



Suf(Su,C, z) dSu dC; 



(9) 



• number of sources brighter than Sn m as a function of 
redshift, 



dN 

dz 



s„>s lh 



J Sum J - c 



f(S„,C,z) dSudC. 



(10) 



3 DATA 

We now describe each of the data sets which we use to con- 
strain the model. 



3.1 Number Counts 

From the choice of possible number counts, we use the fol- 
lowing. 



(i) Spitzer MIPS counts at 70 and 160 |^m (Bethermin 
eTaLpOlO) . 

(ii) BLAST counts from table 3 of P09 (not constrained 
by the CIB) with covariance matrices (upper quadrants of 
tables 4-6 in P09). Since the BLAST counts are given as a 
series of nodes connected by power laws and not counts-in- 
bins, we treat these data differently than counts from other 
instruments. This is discussed further in Section 14.21 We 
have not accounted for the ~ 10 per cent calibration errors, 
which are strongly correlated ( Truch et al.||2009| , or for er- 
rors due to cosmic variance, which may be significant at the 
bright end. 

(iii) Az TEC 1.1mm counts in th e SHADES fields pre- 

( [20101 . A covariance matrix 



sented in Austermann et al. 



is given; however, the correlations are very high, and the 
paper warns about over-interpreting measured correlations, 
so we use only the diagonal elements of the covariance ma- 
trix, ignoring correlations. This may over- weight the AzTEC 
counts, although a simple test (fitting the amplitude of a 
Schechter function with fixed shape) shows that the errors 
are reasonable. 

In Fig. [3] we show the no-evolution counts (Equation [3] 
with the evolutionary parameters p, g and h all set to 1.0). 
These are derived from the 2-parameter local luminosity 
function described in Section [2T] at 250 and 1100 |^m com- 
pared to the measured counts by BLAST and AzTEC, re- 
spectively. We see that while the 250 \im counts are Eu- 
clidean at the bright end, they show strong evolution in the 
Su = 10-100 mjy range. The 1100 |^m counts show strong 
evolution over their full range. 



3.2 Background Intensity 

We use observations of the CIB reported by Fix sen et al.| 
(1998). We choose to fit the model at the same wavelengths 
at which the counts are used: 160, 250, 350, 500, 850 and 
1100 |xm. We assume 30 per cent band-independent errors 



(discussed further in Section 5.3). 



3.3 Redshift Distribution 



We use dN/dz for SCUBA galaxies as measured by Chap- 



man et al. (2005). They present a histogram of 73 sources 



(their fig. 4), but over-plot a model to show 'the likely effects 
of the sample selection.' We assign each of the 9 histogram 
bins a Poisson error, then scale the histogram bin values and 
errors to fit their model. We fit a 3-parameter Gaussian to 
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Figure 3. No-evolution models compared to data at 250 and 
1100 M-m. The error bars at 10~ 4 and 10 1 are part of the BLAST 
250 u.m data set. 



the histogram and use the fitted amplitude to normalise the 
histogram. These scaled bin values and errors are used to 
constrain the shape of the model redshift distribution. We 
assume a limiting flux density Sy im of 5 mJy. 



3.4 Data Sets Not Used 

A number of other relevant data sets exist that, for various 
reasons, we do not use to constrain the model. In most cases, 
the predictions of the best-fit model are compared to the 
unused data sets in Section [5] 



34.1 Counts 

We have chosen to omit the SCUBA 850 \im number counts 
(e.g., |Coppin et al.|2006 ) due to the fact that there is consid- 
erable tension between these measurements and those per- 
formed more recently with the AzTEC camera in a number 
of different fields at 1.1mm. It has been noted since the 
very first surveys were undertaken with AzTEC that in or- 
der to scale its observed galaxy counts to those observed 
with SCUBA, each galaxy, roughly speaking, would have 
to be a factor of ^3 brighter at 850 urn (e.g. |Perera et al.| 



2008), whereas on an object-by-object basis, the individ- 
ual galaxies appear to have a flux density ratio closer to 2, 
which is about what would be expected for a typical SED 
at redshift -2 (e.g. |Perera et al.|2008| |Chapin et"aL|2009 ). 
By itself, this comparison hints at a bias in counting, even 
though individual objects appear to be well-behaved. We 
note that as counting methods have improved, the downward 
correction for biases in the SCUBA counts has increased; 



the Coppin et al. ( 2006 ) measurement falls below essentially 
all of the previous estimates at 5s 50 > 2 mJy, which were 
made using simpler techniques. Similarly, the methodology 
followed by the AzTEC team has also evolved over time and 
more biases have been discovered and corrected tending to 
lower the counts (e.g. |Austermann et al.|2010) Downes et al. 
|2011| . We also note that the counts at 870 [im measured 
with LABOCA fall considerably below those measured with 
SCUBA flWeifi et al.||2009 ). 



We have found that it is impossible to fit a model that 
is consistent with the AzTEC 1.1mm and SCUBA 850 |im 
counts simultaneously within the quoted uncertainties. We 
note that most recent phenomenological models have also 
only attempted to fit SCUBA counts, rather than includ- 
ing other counts near 1 mm (e.g. |Lagache et al.||2003[ |2004| 



Valia nte et al.|2009||Le Borgne et al. 2009; B ethermin et al 
201 1|) . We have explicitly checked that the model of |Valiante 
et al.| ( |2009| , which fits the SCUBA 850 \im counts, signifi- 
cantly exceeds the AzTEC 1.1mm counts, as we find here. 

Taking all of these facts into consideration, there is 
strong evidence that the counts near 1 mm are biased by 
amounts that are not accurately characterised by published 
uncertainties. We feel that the best results will be obtained 
using the more recent, and more sophisticated, AzTEC mea- 
surements. This, however, is clearly an open subject that 
needs to be fully addressed in future assessments. 

Number counts at 1.4 and 2.0 mm measured by the 
South Pole Telescope (SPT) were recently published by 
Vieira et al. ( 2010| ). They measured the bright end of the 
number counts, however, for which contributions by lensed 
galaxies are expected. Since our model does not include lens- 
ing, we do not use these counts to constrain our model. 



3.4-2 Intensity in Redshift Slices 



Pascale et al7] ( |2QQ9[ table 2) presents the CIB at the BLAST 
wavelengths due to 24 \xm sources in six redshift bins from 
z~ 0.4-2. We attempted to use this data set to constrain 
the evolution model by comparing the values to the integral 
of the model over each redshift bin, but had difficulty recon- 
ciling this data set with the others; we found that, particu- 
larly at z < 0.5, the model preferred significantly lower val- 
ues than the BLAST measurements. Since the analysis uses 
photometric redshifts, which may be unreliable, we have not 
included these data in our fits. We also note that the num- 
bers of galaxies in their sample at these lower redshifts are 
extremely small; sampling variance due to clustering is cer- 
tainly a large term that must be added to the quoted Poisson 
uncertainties. 



3.4-3 Redshift Distribution 

A number of measurements of redshifts distributions have 



been made, including at: 170 fxm by ISO (Patris et al.|2003 


Dennefeld et al. 2005 Taylor et al. 


2005); 250, 350 and 


500 \xm by BLAST ( 


Dye et al. 2009 


Dunlop et al. 


2010 


Chapin et al.||2011); and 1.1mm by AzTEC (|Chapin et al. 



2009). We assume limiting flux densities Sn m of 200, 40, 
20 and 3.8 mJy for the 170, 250, 500 and 1100 (xm distribu- 
tions, respectively (we have not considered the 350 \im data 
set). We do not fit these dN/dz distributions, as the selec- 
tion biases in producing the catalogues, due to the need for 
an optical counterpart to identify each source redshift, are 
difficult to quantify. This is also the case for the SCUBA 
distribution that we do use (although an attempt to correct 
the bias has been applied); however, we wanted to include 
at least one data set to constrain the redshift distribution, 
since the model is degenerate without it, and the SCUBA 
sample is the largest available. 
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4 MODEL FITTING 

Given the data sets listed in Section [3] we map out the likeli- 
hood space of the luminosity evolution model using MCMC. 
This allows us to quote most-likely parameter values along 
with errors and correlations between parameters. 



4.1 Parameterisation 

We parameterise the evolution functions p(z) and g(z) as 
connected power laws at a series of nodes at a specific set 
of Zi. We have chosen 6 free parameters in each of p and z 
at a set of z% spaced roughly linearly in log(l +z), spanning 
z = 0-5; the values used are z% =0.1, 0.5, 1.0, 2.0, 3.5 and 5.0. 
The functions p and g are both fixed to 1.0 at z = and to 
10 -12 at z = 7, which serves as a high-redshift cutoff. The 
evolution parameters log(pi) and log (pi) are then linearly 
interpolated in log(l + z). 

4.2 Likelihood Calculations 

We calculate the likelihood of a model with a given set of 
evolution parameters based on % 2 5 



model using 4 chains and run to an 'R-1' convergence of 
0.003 to ensure accurate confidence limits. This typically 
requires 



' 10 5 samples per chain. 



£({&}) = ex P 



1 2 



(11) 



where {£j } represents the N p free parameters and the sum is 
over the data sets. This formulation inherently assumes 
that the errors on the data sets are Gaussian distributed; 
this is not always a valid assumption, but it allows us to 
proceed in a straightforward manner. 

In general, counts are treated as counts-in-bins. At each 
bin centre, we compare the model (Equation |7| to the data. 
The full covariance matrix is used if available. However, as 
previously discussed, the BLAST P(D) counts are treated 
differently: both the model and the data (connected power 
laws) are integrated between each pair of nodes, and the in- 
tegrated values are compared. Errors (including correlations, 
which can be significant) on the data integrals are measured 
using Monte Carlo simulations, sampling from the chains 
produced by P09. We use Gaussians centred on the median 
values and with widths equal to half of the 68 per cent confi- 
dence regions. This is a reasonable description for all values 
except for the integrals between the two faintest bins at each 
wavelength, which have large positive tails. 

For the SCUBA dN/dz, we compare the model to the 
normalised Gaussian discussed in Section 1531 at the six red- 
shift points Zi defined above. Errors at these points are ob- 
tained by propagating errors, including correlations, mea- 
sured by the MCMC. 

We are, in principle, free to set the relative weights w% of 
each data set, but for now the weights are ignored (wi = 1) 
and the relative importance of each data set depends on the 
number of measurements (and their errors) in each set. 

4.3 Monte Carlos 



We use CosmoMCQ as our likelihood sampler (Lewis & Bri- 
dle|2002| ), using Metropolis-Hastings sampling. We run each 



5 RESULTS 

We fit 12 evolutionary parameters {pi and gi at the six zi) 
to the data sets listed in Section [3] We find an important 
dependence on the parameter a, which governs the extent of 



evolution in the colour-luminosity relationship (Section 2.3 ); 
we have tested both a = 1 (colour-luminosity evolution) 
and a = (no colour-luminosity evolution) and find that 
the a = 1 model is a better fit to the data. We therefore 
concentrate on the a = 1 model, but also show results for 
the a — model for comparison. The results of the MCMC 
analysis are presented in Figs.|4| fTo] Throughout, the results 
for the a — 1 model are in blue and for a = in red, with 
the best-fit models in thick solid and thick dashed lines, 
respectively. Realizations of the model, including luminosity 
functions, counts and sample sources lists, are available at 
http : / / cmbr . phas . ubc . ca/ model/ 

The implications of the choice for a, along with fits with 



http : //cosmologist . inf o/cosmomc/ 



a as a free parameter, are discussed further in Section [674] 
Although it is it not shown here, the free a model leads to 
an intermediate value of a = 0.62 zb 0.04, with correlation 
coefficients (to the other parameters) as large as 0.6; the 
correlations are largest for the low-redshift parameters and 
are nearly zero for the high-redshift parameters [zi > 1). 



5.1 Parameters 

The best-fit evolution functions p(z) and g(z), along with 
68 and 95 per cent confidence intervals, are shown in Fig. [4] 
Both models show clear negative density evolution and pos- 
itive luminosity evolution with increasing redshift. Both 
models remain well constrained up to z = 3.5, with the 
exception of the z = 0.5 node in the a — model; this fea- 
ture is discussed further in Section [5.2.11 The data do not 
constrain the models above z — 3.5. 

To guide the eye, we show (as a thin-dashed line) for 
both density and luminosity functions a representative single 
power law, p(z) — (1 + z) 1 and g(z) = (1 + z) s up to a 
cutoff redshift z = 3.5. These values have not been fit to 
the data, but instead fit to the best-fit a = 1 values. The 
power indices of the lines shown are 7 = —6.7 and 5 = 3.5. A 
generic feature of our models is therefore a combined trend of 
negative density evolution and positive luminosity evolution 
as one observes FIR-selected galaxies further into the past. 

Density evolution may be considered an indication of 
the overall galaxy merger rate as a function of time. If galaxy 
formation were to follow a simple bottom-up scenario, such 
as in the case of dark matter halo merger histories, the small- 
est bodies are created first, and over time merge together 
to form a smaller number of galaxies (i.e., positive density 
evolution). Luminosity evolution does not change the total 
number of galaxies in the Universe, but rather their bright- 
ness distribution. Under the previous scenario, one might 
naively expect a large number of less-luminous galaxies to 
merge together in the past, and produce a smaller num- 
ber of more-luminous sources in the present (i.e., negative 
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Figure 4. The parameter values, which are fit to the data sets 
shown in Figs. |6}{To| Above, density evolution, p(z), and below, 
luminosity evolution g(z). The C-L evolution (a = 1) and no C—L 
evolution (a = 0) models are shown, in blue and red, respectively. 
The values of p and g at z = and z = 7 are fixed to 1.0 and 
10 -12 , respectively. The values in between, at 2=0.1, 0.5, 1.0, 
2.0, 3.5 and 5.0, are free parameters in the model. The blue (red) 
symbols connected by solid (long-dashed) lines indicate the best- 
fit parameter values for the a = 1 (a = 0) models. The coloured 
bands represent the 68 and 95 per cent confidence regions. We 
note that pi and gi at each Zi are highly anti-correlated (see 
Fig. [5}. Representative power laws, p(z) = (1 + z) 1 and g(z) = 
(l + z) s up to a cutoff redshift z = 3.5, with 7 = —6.7 and 5 = 3.5, 
have been over-plotted as black short-dashed lines. 



luminosity evolution combined with positive density evolu- 
tion). However, precisely the opposite behaviour is observed 
in deep extra-galactic surveys at different wavelengths. This 



apparent 'cosmic downsizing' (as first noted by Cowie et al. 
1996) has been a topic of great interest to theorists, and our 



modelling results continue to support the trends observed 
at other wavelengths: regardless of the parameterisation, 
there were fewer, but more luminous FIR-bright galaxies 
in the past. For the most luminous of these galaxies to have 
formed in the early Universe, most of their stars must have 
been created in a relatively short burst (less than 1 Gyr for 
galaxies at z i> 3.0, given the available time since the first 
generation of stars during re-ionisation). The FIR- luminous 
galaxies presumably faded slowly over time as the more mas- 
sive, luminous stars completed their life-cycle, giving way to 
lower-mass, cooler and less-luminous stars. 

The parameter distributions and correlations are shown 
in Fig. [5] The distributions for the a = 1 model are shown in 
the lower-left, and for the a = model in the upper-right. 
The diagonal elements show the distributions of the free 
parameters, pi, . . . , p6, #1, • • • , ge- Above and below the di- 
agonals are the joint distributions of all pairs of parameters. 
We note a few features: (i) the zi = 0.5 parameters for the 
a = model are unconstrained, as was noted above; (ii) the 
zq = 5.0 parameters are unconstrained by the data for both 



models, also noted above; (ii) for both p and g, neighbouring 
parameters (e.g., pi and pi+i) are anti-correlated; and (iv) 
pi and gi (at the same Zi) are strongly anti-correlated. 



5.2 Comparison to Data 

Comparisons of the model to the constraining data sets are 
shown in Figs. |6|fTo] 



5.2.1 Number Counts 

Euclidean-normalised differential number counts at a range 
of wavelengths are shown in Figs. [6] and [7] (for the a — 1 
and a — models, respectively). The counts derived from 
the best-fit models are shown as solid lines, with the coloured 
regions showing the 68 and 95 per cent confidence regions. 
Also shown are the contributions to the counts from sources 
binned by redshift, along with the 68 per cent confidence 
regions. The counts at 850, 1400 and 2000 |^m are shown for 
illustrative purposes only, and have not been used in the fits. 

Comparing the counts derived from the a = 1 and 
models, we see that, with the exception of 160 ]im, which the 
a — model is not able to reproduce, the two models give 
very similar counts through the flux density regions covered 
by data. We see, however, that the curves are formed from 
completely different redshift distributions; the a — 1 model 
shows a prominent z — 0.1-0.5 component at all wave- 
lengths, while this component is completely sub-dominant 
in the a — model. This fact provides an explanation for 
why the z = 0.5 density and luminosity evolution param- 
eters in the a — model are completely unconstrained - 
since sources at this redshift do not contribute appreciably 
to the observed counts at any wavelength, there is only an 
upper limit to the strength of sources in this redshift range. 

We remind the reader that the 250, 350 and 500 \im 
BLAST counts are fit using the integral of dN/dS u across 
a set of bins. An example fit (250 \im for the a — 1 model) 
is presented in Fig. [8] This shows very clearly that the two 
models are able to produce essentially the same counts. We 
also point out that the bin-to-bin correlations are in some 
cases quite strong (up to p — 0.95), and that the x 2 calcu- 
lated from the full covariance matrix can be quite different 
from what would be inferred using the diagonal elements 
alone. This is particularly true at 500 |^m, where % 2 =45; if 
we ignore the correlations, we instead find % 2 — 7. 

Although best-fits have been found, it is worth point- 
ing out that none of the overall fits are formally 'good'; in 
particular, the model is quite low at the S u — 10 — 20 mJy 
BLAST counts and is high compared to 1100 \im counts at 
the bright end. We believe that, due to the choice of param- 
eterisation in the P09 P(D) counts, these points are biased 
high. This comes about since a connected power law is not 
able to reproduce the curvature displayed by the model. Re- 
cent results from Herschel/ SPIRE using P(D) to measure 
counts ( Glem Tet al.|2010| , which go to fainter flux densities 
than the BLAST counts, are low compared to BLAST, sup- 
porting this hypothesis (which we discuss further at the end 
of Section 5.3 ). 



Another assumption that may be biasing the fits is 
that the 3 BLAST measurements are correlated with one- 
another, both by instrumental noise and by the fact the same 



10 G. Marsden et al 

-6 -2 -10 10 -2.5 -1.5 -2.1-1.9-1.7-4.0-3.6-3.2 -20 1 3 -10 1.0 1.2 1.4 1.6 1.8 1.9 2.0 2.3 2.5 -20 




-1.5-0.5 0.5 -3 -1 -0.8 -0.4 -2.8-2.4-2.0 -6 -4 -2 -20 0.0 0.4 0.8 0.8 1.2 1.15 1.25 1.5 1.7 2 3 -20 



Figure 5. Likelihood space of log(pi) and log(^) for the two models, as sampled using MCMC. The C-L evolution (a = 1) and no 
C-L evolution (a = 0) models are shown, in blue (lower- left) and red (upper-right), respectively. For each set, the diagonal shows 
the likelihood of each of the free parameters {pi followed by gi, from top-left to bottom-right, for z% = 0.1, 0.5, 1.0, 2.0, 3.5 and 5.0). 
Below/above the diagonal are the correlations between all pairs of parameters. Contours are the 68 and 95 per cent confidence regions. 
We note that p and g at each Zi are highly anti-correlated, and there are moderate anti-correlations between adjacent redshift nodes. 
For both data sets, the z = 5 nodes are not constrained, and the z = 0.5 nodes are unconstrained for the a = model, as can be seen in 
Fig. ^] A spike is apparent in the a = 0, z = 5 nodes; this feature appears to be consistent across and throughout the chains, but as it 
is uncorrelated with the other more well-constrained, lower-redshift parameters, we believe it has no effect on model predictions. 
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Figure 6. Euclidean-normalised differential number counts. The solid lines are the counts derived at each wavelength from the best-fit 
a = 1 model. The coloured bands represent the 68 and 95 per cent confidence regions. The contribution of sources at z = 0.0-0.1, 0.1-0.5, 
0.5-1.5 and 1.5-5.0 are shown as dotted, dashed, dot-dashed and triple-dot-dashed lines, respectively. The coloured bands indicate the 
68 per cent confidence regions. We note that the z = 0.1-0.5 component is sub-dominant at all flux densities and wavelengths, except 
for a narrow region in the shortest bands. The BLAST data sets (250-350 (xm) include poorly-constrained points at low and high flux 
density that are not included in the range shown here; see Fig. |D Counts at 850 urn (SCUBA) and 1.4 and 2.0 mm (SPT) are shown for 
comparison only, and are not used in the fits. 
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Figure 7. Same as Fig. [6] but for a = 0. Note that, with the exception of 160 |U.m, which is not well-fit, the shape of the a = counts 
are very similar to the a = 1 counts (Fig.[6|), but the contribution from the various redshift ranges are completely different. In particular, 
the z = 0.1-0.5 contribution never dominates at any band or flux density. 
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Figure 8. The BLAST 250 u.m P(D) counts integrated between 
nodes compared to models (colour-coded as in the previous fig- 
ures). The horizontal error bars indicate the size of the bins (the 
bin edges correspond to the flux densities of the nodes listed in 
P09 and plotted in Figs. [6] and [7]), while the vertical bars show 
the errors determined from Monte Carlo sampling of the P09 co- 
variance matrix. 



part of the sky has been observed. For both of these reasons, 
it would be desirable to fit the model directly to the maps 
via multi-band P(D). This will be the focus of a later paper. 



5.2.2 Integrated Brightness (CIB) 

The CIB as measured by FIRAS is used as a constraint 
on the integral of intensity over redshift at a range of wave- 
lengths. This is shown in Fig.[9]for both models. The best-fit 
models are shown as filled symbols, with the 68 and 95 per 
cent confidence regions shown as coloured rectangles. We see 
that both models are high compared to FIRAS at nearly all 
wavelengths. It may be interesting to note that the models 



appear to agree with the Lagache et al. ( 2000 ) curve slightly 
better than the Fixsen et al.|(|1998|) curve. 



As with the correlations between bands in the BLAST 
counts, we have also ignored correlations between bands in 
the FIRAS measurement. Proper treatment of these corre- 
lations would likely reduce the overall x 2 of the models. 



5.2.3 Redshift Distribution 



In Fig. 10 we show a variety of dN/dz measurements: 
170 HJn from ISO; 250 and 500 [im from BLAST; 850 [im 
from SCUBA; and 1.1 pun from AzTEC. As discussed in 
Section [3^4] the selection function for these redshift distri- 
butions are very poorly quantified, so we do not make full 
use of these counts as constraints on the model. However, to 
provide at least some direct redshift constraint, we use the 
approximated SCUBA distribution, and at the other wave- 
lengths simply show the counts predicted by the model com- 
pared to the data. 
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Figure 9. The cosmic infrared background as measured by FI- 
RAS ( |Fixsen et al.|[l998| solid line) with representative 30 per 
cent errors (dotted lines). The circular and square points repre- 
sent the CIB derived from the best-fit models, a = 1 and a = 0, 
respectively . The coloured rectangles represent the 68 and 95 per 
cent confidence regions. Also shown, but not used in the fits, are 
the FIRAS (sol id grey line) a nd DIRBE/WHAM measurements 
(diamonds) by Lag ache et al.| ( |2000| . Our models tend to over- 
predict the background values; we discuss the reason for this in 
Section [6] 



5.3 Goodness-of-fit 

Table [l] lists the best-fit % 2 5 along with the contribution 
from each data set, for the a = 1, a = and the free a 
models (columns 3, 9 and 10, respectively). The number of 
data points in each data set is listed in column 2. The total 
X 2 and number of degrees of freedom (DOF) for each fit are 
listed along the bottom two rows. Since there are 62 data 
points and 12 (13) free parameters, there are 50 (49) DOF 
for the fixed (free) a models. The reduced x 2 , X? — xV^dof , 
is an unfortunately high 3.8, 6.8 and 3.5 for the a = 1, a = 
and free a models, respectively. Clearly, the a = 1 model is 
a much better fit to the data than the a = model (but 



see Section 6.5). Adding an extra parameter for the free a 



model also appears to be justified; even so, we focus on the 
a — 1 model, to keep things simple. 

We explore the effects of re-fitting the model omitting 
various data sets in order to probe the 'strain' on the model 
due to any particular data set. We have run the a — 1 model 
on all data sets: (i) excluding 250 and 500 \im counts; (ii) 
excluding 1100 \im counts; (iii) using 850 \im counts instead 
of 1100 \im counts; (iv) excluding the 850 |^m dN/dz distri- 
bution; and (v) excluding the CIB. The results of these tests 
are discussed here: 

(i) The first test (column 4) was meant to explore the ef- 
fects of correlations between the BLAST data sets; however, 
it is not clear how much of the improvement in Xr is due to 
correlations and how much is due to the lessening of tension 
between BLAST and 1100 [im. The value of is slightly 
lower, at 3.1, and we see that the 1100 \im counts and CIB 
are fit much better, although the dN/dz agreement is much 
worse. 

(ii) Removing the 1100 \im counts (column 5) greatly in- 
creases the goodness-of-fit, reducing Xr to 2.2. We see that 
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Table 1. Breakdown of % 2 by data set for different models. The total % 2 and number of degrees of freedom, iVdof, are listed in the 
bottom two rows. In columns 4-8, the a = 1 model is re-fit excluding (or replacing) individual data sets. The x 2 for excluded data sets 
are listed in parentheses. 



BAND No. Pts ol = 1 ol = ol free 

(/xm) All No 250 and No 1100 850 No dN/dz No CIB All All 

data 500 counts counts counts data data 



Counts 



70 


13 


12.2 


14.4 


13.6 


13.5 


14.3 


10.1 


15.6 


14.6 


160 


11 


6.5 


13.9 


7.7 


8.4 


8.9 


5.0 


71.9 


8.3 


250 


6 


10.0 


(28.0) 


13.8 


10.7 


13.1 


11.2 


19.5 


14.0 


350 


5 


16.7 


20.5 


17.3 


17.8 


15.9 


7.7 


47.5 


18.5 


500 


5 


45.3 


(72.7) 


13.8 


25.6 


31.5 


28.5 


81.8 


48.1 


850 


10 


(13.1) 


(22.1) 


(49.1) 


13.7 


(14.6) 


(18.8) 


(20.0) 


(14.8 


1100 


7 


31.6 


7.8 


(916.9) 


(359.0) 


12.0 


26.1 


19.8 


24.7 



Background 



160 


1 


4.0 


2.5 


1.7 


2.5 


1.8 


(3.2xl0 3 ) 


5.1 


1.9 


250 


1 


7.2 


4.0 


4.0 


4.7 


4.3 


(3.9x10 s ) 


7.3 


4.6 


350 


1 


9.6 


3.9 


7.2 


6.9 


6.7 


(2.7xl0 6 ) 


14.2 


7.7 


500 


1 


7.2 


2.0 


7.8 


5.8 


4.6 


(1.7xl0 7 ) 


20.4 


6.4 


850 


1 


1.1 


0.0 


2.9 


1.1 


0.2 


(2.4xl0 7 ) 


11.4 


1.1 


1100 


1 


0.0 


0.4 


0.8 


0.1 


0.1 


(3.6xl0 6 ) 


5.3 


0.1 



850 9 


40.8 


52.8 


dN/dz 
3.3 


7.2 


(238.5) 


24.1 


20.1 


20.8 


Total: 


192.3 


122.2 


93.9 


117.9 


113.4 


112.7 


339.8 


170.6 


N do f- 


50 


39 


43 


53 


41 


44 


50 


49 


X?: 


3.8 


3.1 


2.2 


2.2 


2.8 


2.6 


6.8 


3.5 



the tension between the BLAST and 1100 \im data sets 
is greatly relieved, that the background is reduced at the 
shorter wavelengths, and that dN/dz is allowed to fit nearly 
perfectly. We believe this is a strong clue for developing im- 
proved models, as we discuss in the next section. 

(iii) Fitting the 850 \xm counts instead of the 1100 |^m 
counts (column 6) greatly improves the fit compared to the 
full data set, to a reduced Xr to 2.2. This is because, com- 
pared to 1100 |^m counts, the 850 |^m counts are higher at 
the bright end and lower at the faint end, which allows bet- 
ter fits to the 500 and 850 \im counts. The 850 |xm dN/dz 
distribution is allowed to fit well. We note, however, that the 
model significantly over-predicts the 1.4 and 2.0 mm counts 
(not shown here). 

(iv) Removing the 850 |^m dN/dz constraint (column 7) 
also removes tension, in this case between BLAST and 
1100 urn counts, although not to such a high degree as for 
(ii);here, x ? = 2.8. 

(v) Without the CIB constraints (column 8), we see that 
the integrated background is entirely unbounded. This is 
because the CIB is the only available constraint on the faint 
end of the counts, which dominates the CIB if the faint- 
end counts are sufficiently steep. The quality of fit to all 
other data set is improved, with Xr — 2.6. However, with 
no overall constraint at faint flux densities, the number of 
high-redshift sources is greatly increased. We note that this 
is not reflected in the 850 \im dN/dz constraint, since that 
data set includes only galaxies brighter than Ssso > 5 mJy. 

We have also run a fit of the a = 1 model to test the 
hypothesis that the low-flux density BLAST nodes are bi- 
ased high by the P(D) parameterisation, as discussed in 



Section [5.2. 1[ We have adjusted the 20/15/10mJy nodes of 
the 250/350/500 (xm counts to match the |Glenn et al.| ( |2010 ) 
Herschel/ SPIRE P(D) counts, and have also adjusted the 
values of the other nodes based on the BLAST counts co- 
variance matrices. Additionally, we have doubled the errors 
on the lowest and highest flux density nodes, to compen- 
sate for the non-Gaussian error distributions (faint nodes) 
and cosmic variance (bright nodes). We find that, with this 
adjusted data set, % 2 for the model is 80.3, with Xr — 1-6. 
This model is a much better fit to the 350, 500 and 1100 \im 
counts, the 850 \im redshift distribution, and CIB through 
the BLAST bands (although the model is still ~ 1.5 a high). 
This test indicates that much of the tension in the model is 
due to the P(D) parameterisation, although this alone is not 
the sole cause of the over-predicted CIB. 

To further test how much of this improvement is due to 
the 'adjustment' and how much is due to the doubled error 
bars, we also ran a test with the original BLAST counts, but 
with the errors on the lowest and highest flux density nodes 
doubled. For this test, we find % 2 = 118. This shows that a 
large fraction of the tension between data sets could be due 
to the apparently high data points at low flux density in the 
BLAST P(D) counts, although underestimates of the errors 
at low and high flux density (possibly due to non-Gaussian 
error distributions and cosmic variance, respectively) may 
also be a significant factor. 



6 DISCUSSION 

We now look at some inferences that can be drawn from 
the model, discuss the implications of the colour-luminosity 
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2 3 
Redshift 

Figure 10. The redshift distribution of number counts above a 
flux density limit. The thick solid/long-dashed/dot-dashed curves 
show the predictions of the best-fit a = 1/a = 0/free a models. 
The coloured bands represent the 68 and 95 per cent confidence 
regions. At 170 \im (ISO), 250 and 500 \im (BLAST) and 1.1 mm 
(AzTEC), redshift histograms are shown. At 850 M-m (SCU BA), 
we show the scaled histogram of Chapman et al.l |2b05) with 
Poisson error bars, as described in Section [3. 4. 3| Redshift distri- 
butions at 170 (Patris et al.|2003| [Pennefeld et al.||2005l|Taylo r 



et al.| 2005), 250, 500 M-m jChapin et al.|2011| > and 1.1 mm ( jChapin 



et al. ||2009| > are shown for comparison only, and are not used in 
the fits 



evolution degeneracy and data sets needed to resolve it, and 
consider other possible improvements, including the use of 
new data sets and modifications to the techniques. 



6.1 Colour-colour distributions 



In Section [2 .2| it was shown that our SED library is consis- 
tent with the real spread in galaxy SEDs in the local Uni- 
verse by using them to map our ffiAS-based C) to the 
independent monochromatic luminosity functions at 12, 15, 
25 and 850 \im. However, it is possible that our chosen SED 
shapes could conspire to produce these consistent integral 
quantities, while failing for individual galaxies. Now that we 
have an evolving model of the luminosity function in hand, 
we can go back to our best-fit distribution, apply observa- 
tional selection functions, and compare the distribution of 
colours for individual galaxies from our model to those de- 
tected in real surveys (i.e., to verify the correlations between 
bands). For this test, we have relied on the two best examples 
that we could find in the literature of surveys with colour in- 
formation spanning 24-850 \im: an S70 > 5 mjy Spitzer sur- 




10 100 10 100 

S 70 [mJy] S 60 [Jy] 

Figure 11. Comparison between the colours of galaxies drawn 
from our best-fit evolving luminosity function (contours) with 
data from real surveys (symbols). The left panels are galaxies 
selected at 70 micron above 5mJy (red symbols, |Kartaltepe et ah] 
|2010| . The top-left panel shows 24 \im vs. 70 \im and the bottom- 
left panel shows 160 u.m vs. 70 u.m. There is an approximate flux- 
limit of 36 mJy in the 160 u.m data, which is indicated by a green 
dashed line. The right panels compare our model with SCUBA 
450 and 850 u.m follow-up of 60 |U.m sources brighter than 5.24 Jy 
(blue symbols whose width and height indicate measurement er- 
rors, [Dunn^etl^^ The top right panel 
shows 450 u.m vs. 60 u.m, and the bottom right panel shows 850 u.m 
vs. 60 \im. In all panels, the contours indicate the density of 
sources predicted from the model, with units of (log Jy) -2 deg -2 . 
Unlike Fig. [2] which shows that our SED library and local colour- 
luminosity distribution can produce the correct total numbers of 
galaxies in several bands, this comparison validates our model for 
individual objects. 



vey with cross-matched 24 and 160 |^m observations taken 
as part of the Cosmological Evolution Survey (COSMOS, 
Kartaltepe et al. 2010); and the Submillimetre Local Uni- 
verse Galaxy Survey (SLUGS), in which SCUBA was used 
to follow up #60 > 5.24 Jy IRAS galaxies at 450 and 850 |^m 
( |Dunne et al.|2000||Dunne fe Eales|2001) . Almost all of the 



SLUGS galaxies have luminosities L > 1O 1O L , with about 
half above L > 10 11 L©, and it is a truly local sample, with 
all of the objects lying at z < 0.1. The COSMOS sample is 
deeper and higher-redshift, though with a similar range of 
luminosities; about 50% of the sources lie at z > 0.5, and 
virtually all of the objects have luminosities L > 1O 1O L0, 
while rsj 70% have luminosities L > 10 11 L©. 

The comparison of our model to these data sets is 
shown in Fig. [TT] The Spitzer COSMOS data shows ex- 
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cellent correspondence with our model in terms of the 160- 
to-70 \xm colour. Our theoretical distribution also broadly 
reproduces the main trend in the observed 24-to-70 \im cor- 
relation, although the scatter is slightly less than that ob- 
served in real galaxies. As stated in the introduction, we 
have not attempted to fully reproduce the observed prop- 
erties of 24 urn sources, since the SEDs of galaxies at those 
shorter wavelengths depend on other factors (which would 
then require more variables). We also note that since this 
is a higher-redshift sample, the 24 urn band will sample the 
PAH features for a small subset. The correspondence with 
the SLUGS galaxies is excellent, showing that the spread 
and correlation between the Ray leigh- Jeans side of our SED 
templates and the FIR peak is accurate. 



6.2 Star Formation Rate History 

Star formation rate density (SFRD) can be derived from the 
model by calculating the luminosity density Ltik(z): 



dL 



dC 



L d ^{C) 



9(L,C,z), (12) 



where Ltir is the luminosity integrated over the total IR 
range, 8-1000 |^m, and dLTm/dL, which depends only on 
C, converts L to Ltir- This can then be converted to SFRD 
using the simple relation given by |Kennicutt] ( |1998| , 

-10, 



SFR [M yr" 1 ] = 1.728 x 10" 1u L T ir [L ] 



(13) 



The resulting SFRD for the two models are shown in Fig. |12| 
The shapes are very different, with the a = 1 model peaking 
at z = 1 and the a = peaking at z = 2, but with no SFRD 
whatsoever at z < 1. We believe that this is due to the data 
requiring a certain amount of cool- type galaxies; this can 
either come from intrinsically cool SEDs or from redshifting 
moderate-temperature SEDs so that they appear cooler. The 
a — 1 model provides a higher fraction of cool galaxies with 
increasing redshift, while the a — model is required to 
compensate by strongly increasing the number of galaxies 
at high redshift; it is then necessary to decrease the number 
of galaxies at low-to-moderate redshift. The free a model is 
also shown (in green, dash-dotted line). Since the best-fit a 
value is ^0.6, it is not surprising that the SFRD curve falls 
between the a — 1 and models. 

We show a small sampling of data along with the mod- 
els. Below z ~ 2, the data are generally consistent with each 
other and are broadly in line with the a — 1 model; this 
model in particular agrees well with the Wall et al. (2008) 
compilation. Above z ~ 2, the data are inconsistent and it is 
hard to draw any firm conclusions, although we note that the 
|Steidel et al.| ( |1999| and |Giavalisco et al.| ( |2004| points are 
based on extinction-corrected UV measurements and do not 
necessarily bear any relation to SFRD estimates in the FIR. 
The a = model appears to be inconsistent with the data. 
We also note that the model with the 'adjusted' BLAST 
counts, as described at the end of Section |5.3| shows an 
SFRD curve that is slightly higher at z = 0.5, slightly lower 
at z = 1, and runs a bit flatter between z = 2 and 3.5 com- 
pared to the displayed a = 1 model, fitting the |Wall et ah] 
(2008) compilation slightly better. 




2 3 4 

Redshift 

Figure 12. Star formation rate density as a function of redshift. 
The a = 1 (blue), a = (red) and free— a (green) models are 
shown. A sample of other estimated data points are shown; data 



are from|Pascale et al. (2009), Wall et al. 


2008), Giavalisco et al. 


(2004 


|,|Lilly et al.| (| 1 996 |).| Connolly et al. 


1997J and|Steidel et al. 


(1999 


I (the latter four taken fromlMichalowski et al.|2010| which 



has an extensive list of values in its appendix). Note that the 
jagged appearance in the model predictions is partly a result of 
the anti-correlations between adjacent redshift nodes in the fits 
(as seen in Fig. [5j. 



6.3 Fraction of CIB Resolved By BLAST 

If the CIB is in fact as high as indicated by the models (^2a 
high compared to FIRAS at the BLAST wavelengths for the 
a = 1 model, Fig. |9j), the fraction of the CIB resolved by 
stacking the BLAST maps on 24 \im Spitzer sources would 
in fact be lower than is quoted by Marsden et al. (2009). 



Assuming the CIB values predicted by the a — 1 model, 
I u — 1.5, 1.3 and 0.7MJysr _1 , the stacked intensities mea- 
sured by Marsden et al.| ( |2009| correspond to 45, 45 and 50 
per cent of the CIB, respectively. We note, however, that 
the CIB predicted by the model is likely driven high by the 
apparently high BLAST nodes, discussed in Section |5.2.1| 
Comparing to the model with the modified BLAST points, 
we find 55, 55 and 65 per cent, still significantly lower than 



the percentages quoted in Marsden et al. (2009), due to the 
fact that the model CIB is still ~ 1.5 a high through the 
BLAST bands. 



6.4 Colour-Luminosity Evolution 

We have shown that, based on counts alone, the a = 1 and 
a — models are nearly indistinguishable. The redshift dis- 
tributions, however, are very different, with the a — model 
a better fit to the scaled [Chapman et al. ( 2005| 850 \xxn. red- 
shift distribution. The SFRD figure, in particular, shows a 
large discrepancy between the models, with the data signif- 
icantly preferring the a = 1 model. To truly rule out one 
model or the other, however, direct measurements of the 
redshift distribution are needed. 

In Fig. [13] we show the integrated 350 \im counts of 
sources brighter than 20 mjy in the colour-redshift plane for 
both models, with one-dimensional distributions shown in 
the side panels. The low-redshift distributions are of course 
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Figure 13. The density of 350 \±m sources brighter than 20mJy 
in the colour-redshift plane for both a = 1 (blue, solid contours) 
and a = (red, long-dashed contours) models. The background 
intensity maps are linear and the contours are 1, 10, 100 and 
1000 deg -2 per unit colour per unit redshift. The distributions 
marginalised over redshift and colour are shown at the right and 
top, respectively. 



exactly the same, but, as previously seen, the models are 
distinct at higher redshifts. We see that the a = 1 model 
peaks at C ~ — 0.4, z~ 1, while the a = model peaks at 
C ~ 0.1, z~2. A SPIRE redshift survey down to a flux limit 
of 20mJy will allow us to unambiguously disentangle the 
two models. 

When attempting to interpret the results of submm 
flux-limited samples, we note the strong bias toward detect- 
ing cooler, less-luminous sources. Following the prescription 



of Chapin et al. (2011 ), we calculate the probability density 
in the colour- luminosity plane at two fixed redshifts, z = 0.1 
and 1.0, for samples selected at 350 \im above a flux density 
of 20 mJy, for both the a = 1 and a = models. The results 
are shown in Fig. |14[ For reference, the locus of the 350 |im 
flux limit is shown as short- dashed lines at each redshift, and 
the thin solid line shows the local colour-luminosity correla- 
tion. The sample is biased downwards and to the left; this 
is due to the fact that the total number of sources increases 
toward lower luminosities, and that the selection function is 
inclined in the colour-luminosity plane; cooler sources may 
be detected at fainter luminosities, where there are consider- 



ably more objects. IChapin et al. (2011) used this argument 



to suggest that their submm-selected sample could be consis- 
tent with non-evolution of the colour-luminosity correlation 
(despite finding a tail of cool sources below the local colour- 
luminosity correlation). Similarly, included among the Sci- 
ence Demonstration Phase results from Herschel are some 
tentative indications that submm-selected sources are cooler 
in the past (e.g., Hwang et al. 2010 Elbaz et al. 2010); it 



will be necessary to model this selection effect to determine 
whether these observed trends are indeed real. 



Of course, another possibility is that the colour- 
luminosity correlation evolves in a way completely different 
from the simple form assumed in this paper. For example, 
Symeonidis et ah] ( |2010| argue that the correlation simply 
broadens at high redshift, based on a sample of galaxies 
spanning the submm-FIR, rather than there being a trend 
toward lower temperatures in the past. While the observed 
correlation at high-redshift may in fact be broader than the 
local correlation, we note that the combination of a tight 
colour-luminosity correlation with a submm selection effect 
can also broaden the apparent range of colours, as shown in 
Fig. [13] and certainly contributes, at least in part, to the 



effect seen by Symeonidis et al.|(2010) 



Based on the forward-modelling approach taken in this 
paper, which intrinsically accounts for these selection bi- 
ases, we find that an evolution toward cooler temperatures 
in the past is the most likely scenario (i.e., the tail of cooler 
sources in the high-redshift submm-selected population is 
even cooler than what one might expect given the selec- 
tion biases) . This conclusion is in rough agreement with the 
model of Valiante et al. (2009}, who also included luminos- 



ity evolution in the colour-luminosity correlation to produce 
cooler SEDs in the past, although they also included an 
extra population of cold local sources that are presumably 



missed in IRAS surveys. In contrast, Lewis et al. (2005) 



and |Le Borgne et al. ( 2009 ) used fixed colour-luminosity 
correlations as a function of redshift, and also appear to 
obtain reasonable results. We note, however, that our re- 
sult is based strongly on the imposed colour distribution at 
the faint end of the luminosity function, which we believe 
is possibly the cause of the over-predicted CIB. A modified 
colour-luminosity relationship could have a large effect on 
the relative merits of the a — 1 and a — models. This is 
discussed further in the next section. 

As already discussed, the nature of the evolution of the 
colour-luminosity correlation has a direct impact on our in- 
ferences about the total SFRD as a function of time. As 
shown in Fig. [12] a decrease in the typical temperature of 
galaxies at high redshift leads to a later peak formation 
epoch when compared with a non-evolving colour-luminosity 
scenario. Working from primarily SCUBA-selected samples, 
there is some evidence that FIR-luminous galaxies at high- 
redshift are indeed cooler, but also physically more ex- 



tended, based on radio morphologie s (e.g., |Chapman et al 
|2004|), near-mid-IR colours (e.g., |Hainline et al.||2009|), 



and mid-IR spectra (e.g., Men endez-Delmestre et al.|2009 ). 
These observations appear to be consistent with a local- 
Universe measurement that showed an anti-correlation be- 
tween physical size as a function of luminosity and temper- 
ature ( |Chanial et al.||2007[ ). 



6.5 Reconciling the CIB with the SMG z > 1 
redshift distribution 

As has been seen in previous sections, the single-population 
model that we have attempted to fit fails, primarily, in rec- 
onciling the redshift distribution of ^1 mm-selected SMGs 
(and therefore the z > 1 SFRD, Figs. 



10 



and 



12), with the 



intensity of the CIB (Fig. [9]), despite doing a decent job of fit- 
ting the counts in the bands that we have considered (Figs. [6] 
and|7|. Generally speaking, we have found that models with 
lower (and therefore more consistent) values of the CIB push 
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frared galaxies (L > 1O 12 L0) are merger- driven starbursts 
( Sanders & Mirabel 1996 ). In contrast, the faint end must in- 



Luminosity [L©] 

Figure 14. The density of 350 u.m sources brighter than 20 mJy 
at a particular redshift in the colour-luminosity plane for both 
a = 1 (blue, solid contours) and a = (red, long-dashed con- 
tours) models. We show the population at two redshifts, z = 0.1 
and 1.0. The contours correspond to 0.1 and 0.5 of the maxi- 
mum density. The thin short-dashed lines show the flux limits at 
the two redshifts. These lines are independent of luminosity be- 
low C = —0.5 due to our choice of SED library. The solid black 
line shows the local colour-luminosity relationship, p(C\L) (see 
Section [2.1| , with 1-cr limits plotted as dotted lines. This figure 
demonstrates the strong bias toward cooler and lower-luminosity 
sources when compared to the rest-frame distribution, caused by 
our selection function. 



the ^lmm selected galaxies to lower (and less consistent) 
redshifts. 

We believe most of the discrepancy is due to our lack 
of knowledge of the low- luminosity (L < 1O 1O L0) galaxy 
SEDs. These faint galaxies are essentially undetected indi- 
vidually above the confusion limit in any of the existing 
blank- field ground-based ~1 mm surveys (Blai n et al.|2002|), 
nor in any BLAST 250-500 \im (|Dye et al.||2009| |ChaphT 



et al. 2011), or Spitzer 70 and 160 jam surveys (see, for 



example, the description of the COSMOS survey in Sec- 
tion 6.1 ). Even IRAS struggled to conduct unbiased surveys 



of such faint objects, since it was limited in sensitivity and 
also the nearby galaxies in the local over-density were lim- 



ited by sampling variance (see Section 2.1). However, the 
integrated light from these faint but numerous galaxies can 
make a large contribution to the CIB. 

Referring to Fig. 4 of C09a, it is clear that the clean 
correlation between C and L broadens significantly between 
10 10 L to 1O 9 L , and then it completely breaks down at 
even lower luminosities. Much of this effect is probably due 
to the fact that the 60/100 ]im flux ratio is simply a poor 
indicator of the longer-wavelength SED of galaxies with ex- 
tremely cool dust (i.e., these bands sample relatively warmer 
and potentially un-related dust). It is probably for this rea- 
son that other authors have found it necessary to modify 
by hand the distribution of low-luminosity / cool galaxies in 
their evolutionary models. 

In addition, the faint end of the FIR luminosity func- 
tion contains a significantly more heterogeneous collection 
of galaxies than the bright end. Most luminous infrared 
galaxies (L > 1O 11 L ) and virtually all ultra-luminous in- 



clude everything else, including mildly star-forming spirals 
like the Milky Way (relatively low-luminosity despite the 
presence of dust), passively-evolving spheroids (extremely 
faint due to the near complete lack of dust), but also smaller 
star-forming galaxies; blue compact dwarfs, for example, 
have dust and contain a number of active star-forming re- 
gions which may produce apparently large dust tempera- 
tures despite their lower luminosities. 

To test the hypothesis that our lack of knowledge at 
these low luminosities plays a crucial role, we performed a 
simple test using our best-fit a = 1 model. First, we assigned 
all galaxies below 10 10 L the warmest SED (largest value 
of C). This has the effect of decreasing the flux densities at 
A > 100 |^m, and increasing the flux densities at A < 60 |^m 
when extrapolating from the FIR luminosity for those galax- 
ies. We then re-calculated the number counts, CIB spec- 
trum, and redshift distribution of 850 |^m selected galaxies. 
While the modification to the SEDs of these faint galaxies 
had little impact on either the counts or redshift distribu- 
tion, where data are available (again, since galaxies with 
such luminosities lie below current confusion limits), there 
was a huge change to the CIB. Its predicted spectrum be- 
came much warmer (peaking closer to 100 fim than 200 |im), 
and dropped significantly below the data at submm wave- 
lengths. Conversely, assigning these low-luminosity galaxies 
the coolest SED has the opposite effect, pushing the pre- 
dicted CIB peak to longer wavelengths, and exceeding the 
CIB by even more than with the regular model. 

We then repeated this procedure using a lower luminos- 
ity threshold of 10 9 L . While the sense of the changes to the 
predicted quantities were the same, the impact is obviously 
smaller, and using the warmest SEDs for these faint galax- 
ies, we obtain a reasonable fit to the CIB (slightly low at 
200 \im, and then going through the data at longer submm 
wavelengths). 

Based on this experiment, it is clear that how one treats 
galaxies at the faint-end of the FIR luminosity function is 
crucial, even although they are not well constrained observa- 
tionally. While one could modify a phenomenological model 
ad hoc to evolve the faint-end significantly less than the 
bright-end as in Lagache et al. (2003) and Valiante et al. 
( 2009 ) , thus reducing the number of cool-galaxies at higher 
redshifts, the alternative of increasing the temperatures of 
fainter galaxies at all redshifts would have a similar effect. In 
the future one could use wide-area surveys such as H- ATLAS 



Eales et al. (2010) to measure the SEDs of fainter nearby 
galaxies, and perhaps also with SCUBA-2 at 450 p,m(since 
it will be capable of resolving most of the CIB directly into 



individual sources; Holland et al. (2006)). We also note that 
since our present model tends to pull down the SMG redshift 
distribution in order to improve the \ 2 of the CIB, our re- 
sult that the a = 1 model is better than the higher-redshift 
a = model is far from secure. Soon it should be able to 
perform tests such as Figs. [!3}fT4] to determine, at least for 
the more luminous objects, whether there really is evolu- 
tion in the luminosity-temperature correlation using SPIRE 
data, particularly from HerMES. 
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6.6 Future Improvements 

Although in many ways our model is an improvement over 
previously published studies, throughout the paper we have 
discussed several shortcomings of the model and data sets 
used. Here we list some analysis techniques and future data 
sets that will improve the quality of the model: 



• As discussed in Section |5TT] the BLAST P(D) counts 
are inherently correlated. Instead of using the counts as an 
intermediate data set, we can use P{D) within the model- 
fitting framework to fit the maps directly. This is of course 
computationally intensive, but certainly worth pursuing. It 
also has the advantage over the P09-style P(D) in that the 
shape of the model is much more reflective of the shape of the 
true underlying counts, rather than imposing a connected 
power law onto the counts. Furthermore, as we noted in 
Section [5.4.1 1 the discrepancy between SCUBA 850 \im and 
AzTEC 1.1 mm counts may be due to the potentially biased 
counting of sources. P{D) analyses do not suffer most of 
those problems. 

• New redshift distributions from SPIRE will go a long 
way to constraining the model, in particular breaking the 
colour-luminosity evolution degeneracy, and enabling more 
redshift nodes to be used. 

• Direct measurements of the local luminosity functions 
at FIR and submm wavelengths will give us a much better 
starting point for the evolving luminosity function, and will 
allow us to derive a better SED library. Wide-area PACS 
and SPIRE surveys will be able to provide these. We also 
hope that wide-area surveys will be capable of detecting 
and measuring the SEDs of galaxies with FIR luminosities 
in the range 10 10 -10 11 to solve the discrepancy between 
the redshift distribution of SMGs and the spectrum of the 
CIB. 

• Extend the model to the mid-IR by including a more 
sophisticated SED library. This may require additional SED 
parameters to include, for example, the AGN contribution 
(e.g. Valiante et al.|2009| ). 

• Include lensing by adopting a similar approach to 
Paciga et al.| ( |2009| . Such a treatment is now necessary to 



explain the counts at mm wavelengths covering wide areas, 
and more recently, Herschel/ SPIRE surveys. 

• More versatile modelling of the evolving SED distribu- 
tion, in particular changing the width and shape of P(C\L) 
as a function of L. 

Finally, we note that a major part of the work in this pa- 
per went into developing likelihood expressions for the var- 
ious data sets. To improve on our methodology for current 
and future surveys, it will be important to fully characterise 
uncertainties and correlations between data sets. For exam- 
ple, the differential number counts, integrated background 
and redshift distributions of sources in the same field will 
have a correlated cosmic variance term. 



7 CONCLUSIONS 

We have presented a sophisticated technique using MCMC 
to fit a simple evolving luminosity function to a range of 
FIR and submm data. We are able to measure errors on the 
evolutionary parameters and show the correlations between 



these parameters. The results of the model are available at 
http : //cmbr . phas . ubc . ca/model/ We also show how the 
various data sets are in tension with one another and demon- 
strate the importance of redshift distributions. 

An advantage of our approach over some other models 
in the literature is that we need only consider the evolution 
of one galaxy population with a single-parameter family of 
SEDs based on the correlation between the 60-to-100 \im 
rest-frame colour and FIR luminosity. While we find that, 
across the 70-1100 |^m wavelength range, the counts can be 
fit using models with and without evolution toward cooler 
galaxy dust temperatures at higher-redshifts, there is sig- 
nificant tension between the spectrum of the CIB and the 
redshift distribution of SMGs. We believe that most of this 
discrepancy is caused by the presently unknown distribution 
of submm SEDs for galaxies with luminosities < 1O 1O L0. 
Emerging data from Herschel will immediately help in two 
areas: (i) wide- area surveys, such as H- ATLAS, should en- 
able us to measure the SEDs for at least a small sample 
of such objects; (ii) deeper surveys, such as HerMES, can 
potentially be used to search for evolution in the colour- 
luminosity correlation at the brighter end of the luminosity 
function, as indicated in Figs. |13f|14| Such data should obvi- 
ate the need for ad hoc modifications to the low-luminosity 
region of the local luminosity function by, for example, us- 
ing multiple uncorrelated galaxy populations. The reality is 
that there is a continuum of galaxy types, and it is our hope 
that we can gain more realistic insight into how they form 
and evolve using the simplest phenomenological models that 
are consistent with the data. 
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